function [num, den, z, p] = besselfilter(n, Wn, ftype, anaflag)
%BESSELF  Bessel analog filter design.
%   [B,A] = BESSELF(N,Wn) designs an N'th order lowpass analog
%   Bessel filter and returns the filter coefficients in length
%   N+1 vectors B and A.  The cut-off frequency Wn must be
%   greater than 0.
%
%   If Wn is a two-element vector, Wn = [W1 W2], BESSELF returns an 
%   order 2N bandpass filter with passband  W1 < W < W2.
%   [B,A] = BESSELF(N,Wn,'high') designs a highpass filter.
%   [B,A] = BESSELF(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2].
%   
%   When used with three left-hand arguments, as in
%   [Z,P,K] = BESSELF(...), the zeros and poles are returned in
%   length N column vectors Z and P, and the gain in scalar K. 
%
%   When used with four left-hand arguments, as in
%   [A,B,C,D] = BESSELF(...), state-space matrices are returned.
%
%   See also BESSELAP, BUTTER, CHEBY1, CHEBY2, FREQZ and FILTER.

%   References:
%     [1] T. W. Parks and C. S. Burrus, Digital Filter Design,
%         John Wiley & Sons, 1987, chapter 7, section 7.3.3.

btype = 1;
if (nargin >= 3)
% band-stop or high-pass
	btype = 3; 
end
if length(Wn) == 2
	btype = btype + 1;
end

% step 1: get analog, pre-warped frequencies
u = Wn;

% step 2: convert to low-pass prototype estimate
if btype == 1	% lowpass
	Wn = u;
elseif btype == 2	% bandpass
	Bw = u(2) - u(1);
	Wn = sqrt(u(1)*u(2));	% center frequency
elseif btype == 3	% highpass
	Wn = u;
elseif btype == 4	% bandstop
	Bw = u(2) - u(1);
	Wn = sqrt(u(1)*u(2));	% center frequency
end

% step 3: Get N-th order Butterworth analog lowpass prototype
[z,p,k] = besselap(n);

% Transform to state-space
[a,b,c,d] = zp2ss(z,p,k);

% step 4: Transform to lowpass, bandpass, highpass, or bandstop of desired Wn
if btype == 1		% Lowpass
	[a,b,c,d] = lp2lp(a,b,c,d,Wn);

%elseif btype == 2	% Bandpass
%	[a,b,c,d] = lp2bp(a,b,c,d,Wn,Bw);

%elseif btype == 3	% Highpass
%	[a,b,c,d] = lp2hp(a,b,c,d,Wn);

%elseif btype == 4	% Bandstop
%	[a,b,c,d] = lp2bs(a,b,c,d,Wn,Bw);
end

if nargout == 4
	num = a;
	den = b;
	z = c;
	p = d;
else	% nargout <= 3
% Transform to zero-pole-gain and polynomial forms:
	if nargout == 3
		[z,p,k] = ss2zp(a,b,c,d,1);
		num = z;
		den = p;
		z = k;
	else % nargout <= 2
		den = poly(a);
                [z,k] = tzero(a,b,c,d);
                num = k * poly(z);
                num = [zeros(1,length(den)-length(num))  num];
	end
end


function [a,b,c,d] = zp2ss(z,p,k)
%ZP2SS	Zero-pole to state-space conversion.
%	[A,B,C,D] = ZP2SS(Z,P,K)  calculates a state-space representation:
%		.
%		x = Ax + Bu
%		y = Cx + Du
%
%	for a system given a set of pole locations in column vector P,
%	a matrix Z with the zero locations in as many columns as there are
%	outputs, and the gains for each numerator transfer function in
%	vector K.  The A,B,C,D matrices are returned in block diagonal
%	form.
%
%   	This function handles SIMO systems if the Control Toolbox is
%   	present and SISO systems if only the Signal Processing Toolbox
%   	is installed.
%   	
%	See also SS2ZP.


[mn,nn] = size(z); [md,nd]=size(p);
% Put it in column format if its SISO and in row format.
if (length(k)==1 & md < 2 & mn < 2) & (nn > 1 | nd > 1) 
	z = z'; p = p'; 
end
[m,n] = size(z);
if n==0, n=length(k); end  % Fix to handle multi-output when z is empty
if length(k) ~= n & (~isempty(z))
	error('Z should be a column vector or K should be SIMO.');
end
if n > 1
	% If it's multi-output, we can't use the nice algorithm
	% that follows, so use the numerically unreliable method
	% of going through polynomial form, and then return.
	eval('[num,den] = zp2tf(z,p,k);') % Suppress compile-time diagnostics
	[a,b,c,d] = tf2ss(num,den);
	return
end

% Strip infinities and throw away.
p = p(finite(p));
z = z(finite(z));

% Group into complex pairs
np = length(p);
nz = length(z);
z = cplxpair(z,1e6*nz*norm(z)*eps + eps);
p = cplxpair(p,1e6*np*norm(p)*eps + eps);

% Initialize state-space matrices for running series
a=[]; b=[]; c=[]; d=1;

% If odd number of poles AND zeros, convert the pole and zero
% at the end into state-space.
%	H(s) = (s-z1)/(s-p1) = (s + num(2)) / (s + den(2))
if rem(np,2) & rem(nz,2)
	a = p(np);
	b = 1;
	c = p(np) - z(nz);
	d = 1;
	np = np - 1;
	nz = nz - 1;
end

% If odd number of poles only, convert the pole at the
% end into state-space.
%  H(s) = 1/(s-p1) = 1/(s + den(2)) 
if rem(np,2)
	a = p(np);
	b = 1;
	c = 1;
	d = 0;
	np = np - 1;
end	

% If odd number of zeros only, convert the zero at the
% end, along with a pole-pair into state-space.
%   H(s) = (s+num(2))/(s^2+den(2)s+den(3)) 
if rem(nz,2)
	num = real(poly(z(nz)));
	den = real(poly(p(np-1:np)));
	wn = sqrt(prod(abs(p(np-1:np))));
	if wn == 0, wn = 1; end
	t = diag([1 1/wn]);	% Balancing transformation
	a = t\[-den(2) -den(3); 1 0]*t;
	b = t\[1; 0];
	c = [1 num(2)]*t;
	d = 0;
	nz = nz - 1;
	np = np - 2;
end

% Now we have an even number of poles and zeros, although not 
% necessarily the same number - there may be more poles.
%   H(s) = (s^2+num(2)s+num(3))/(s^2+den(2)s+den(3))
% Loop thru rest of pairs, connecting in series to build the model.
i = 1;
while i < nz
	index = i:i+1;
	num = real(poly(z(index)));
	den = real(poly(p(index)));
	wn = sqrt(prod(abs(p(index))));
	if wn == 0, wn = 1; end
	t = diag([1 1/wn]);	% Balancing transformation
	a1 = t\[-den(2) -den(3); 1 0]*t;
	b1 = t\[1; 0];
	c1 = [num(2)-den(2) num(3)-den(3)]*t;
	d1 = 1;
%	[a,b,c,d] = series(a,b,c,d,a1,b1,c1,d1); 
% Next lines perform series connection 
	[ma1,na1] = size(a);
	[ma2,na2] = size(a1);
	a = [a zeros(ma1,na2); b1*c a1];
	b = [b; b1*d];
	c = [d1*c c1];
	d = d1*d;

	i = i + 2;
end

% Take care of any left over unmatched pole pairs.
%   H(s) = 1/(s^2+den(2)s+den(3))
while i < np
	den = real(poly(p(i:i+1)));
	wn = sqrt(prod(abs(p(i:i+1))));
	if wn == 0, wn = 1; end
	t = diag([1 1/wn]);	% Balancing transformation
	a1 = t\[-den(2) -den(3); 1 0]*t;
	b1 = t\[1; 0];
	c1 = [0 1]*t;
	d1 = 0;
%	[a,b,c,d] = series(a,b,c,d,a1,b1,c1,d1);
% Next lines perform series connection 
	[ma1,na1] = size(a);
	[ma2,na2] = size(a1);
	a = [a zeros(ma1,na2); b1*c a1];
	b = [b; b1*d];
	c = [d1*c c1];
	d = d1*d;

	i = i + 2;
end

% Apply gain k:
c = c*k;
d = d*k;

return

function [a,b,c,d] = tf2ss(num, den)
%TF2SS	Transfer function to state-space conversion.
%	[A,B,C,D] = TF2SS(NUM,DEN)  calculates the state-space 
%	representation:
%		.
%		x = Ax + Bu
%		y = Cx + Du
%
%	of the system:
%		        NUM(s) 
%		H(s) = --------
%		        DEN(s)
%
%	from a single input.  Vector DEN must contain the coefficients of
%	the denominator in descending powers of s.  Matrix NUM must 
%	contain the numerator coefficients with as many rows as there are
%	outputs y.  The A,B,C,D matrices are returned in controller 
%	canonical form.  This calculation also works for discrete systems.
%	To avoid confusion when using this function with discrete systems,
%	always use a numerator polynomial that has been padded with zeros
%	to make it the same length as the denominator.  See the User's
%	guide for more details.
%
%	See also: SS2TF.

[mnum,nnum] = size(num);
[mden,n] = size(den);
% Check for null systems
if  (n == 0 & nnum == 0), a=[]; b=[]; c=[]; d=[]; return, end

% Strip leading zeros from denominator
inz = find(den ~= 0);
den = den(inz(1):n);
[mden,n] = size(den);

% Check for proper numerator
if nnum > n
	% Try to strip leading zeros to make proper
	if (all(all(num(:,1:(nnum-n)) == 0)))
		num = num(:,(nnum-n+1):nnum);
		[mnum,nnum] = size(num);
	else
		error('Denominator must be higher or equal order than numerator.');
	end
end

% Pad numerator with leading zeros, to make it have the same number of
% columns as the denominator, and normalize it to den(1)
num = [zeros(mnum,n-nnum) num]./den(1);

% Do the D-matrix first
if length(num)
	d = num(:,1);
else
	d = [];
end

% Handle special constant case:
if n == 1
	a = [];
	b = [];
	c = [];
	return
end

% Now do the rest, starting by normalizing den to den(1),
den = den(2:n) ./ den(1);
a = [-den; eye(n-2,n-1)];
b = eye(n-1,1);
if mnum > 0
	c = num(:,2:n) - num(:,1) * den;
else
	c = [];
end

function [num,den] = zp2tf(z,p,k)
%ZP2TF	Zero-pole to transfer function conversion.
%	[NUM,DEN] = ZP2TF(Z,P,K)  forms the transfer function:
%
%	                NUM(s)
%	        H(s) = -------- 
%	                DEN(s)
%
%	given a set of zero locations in vector Z, a set of pole locations
%	in vector P, and a gain in scalar K.  Vectors NUM and DEN are 
%	returned with numerator and denominator coefficients in descending
%	powers of s.  
%

%	Note: the following will not work if p or z have elements not
%	in complex pairs.

den = real(poly(p(:)));
[md,nd] = size(den);
k = k(:);
[mk,nk] = size(k);
if isempty(z), num = [zeros(mk,nd-1),k]; return; end
[m,n] = size(z);
if mk ~= n
	if m == 1
		error('Z and P must be column vectors.');
	end
	error('K must have as many elements as Z has columns.');
end
for j=1:n
	zj = z(:,j);
	pj = real(poly(zj)*k(j));
	num(j,:) = [zeros(1,nd-length(pj)) pj];
end

function [zd, pd, kd, dd] = bilinear(z, p, k, fs, fp, fp1)
%BILINEAR Bilinear transformation with optional frequency prewarping.
%   [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs) converts the s-domain transfer
%   function specified by Z, P, and K to a z-transform discrete
%   equivalent obtained from the bilinear transformation:
%
%      H(z) = H(s) |
%                  | s = 2*Fs*(z-1)/(z+1)
%
%   where column vectors Z and P specify the zeros and poles, scalar
%   K specifies the gain, and Fs is the sample frequency in Hz.
%   [NUMd,DENd] = BILINEAR(NUM,DEN,Fs), where NUM and DEN are 
%   row vectors containing numerator and denominator transfer
%   function coefficients, NUM(s)/DEN(s), in descending powers of
%   s, transforms to z-transform coefficients NUMd(z)/DENd(z).
%   [Ad,Bd,Cd,Dd] = BILINEAR(A,B,C,D,Fs) is a state-space version.
%   Each of the above three forms of BILINEAR accepts an optional
%   additional input argument that specifies prewarping. For example,
%   [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs,Fp) applies prewarping before
%   the bilinear transformation so that the frequency responses
%   before and after mapping match exactly at frequency point Fp
%   (match point Fp is specified in Hz).
%
%   See also IMPINVAR.

%   Gene Franklin, Stanford Univ., motivated the state-space
%   approach to the bilinear transformation.

[mn,nn] = size(z);
[md,nd] = size(p);

if (nd == 1 & nn < 2) & nargout ~= 4	% In zero-pole-gain form
	if mn > md
		error('Numerator cannot be higher order than denominator.')
	end
	if nargin == 5		% Prewarp
		fp = 2*pi*fp;
		fs = fp/tan(fp/fs/2);
	else
		fs = 2*fs;
	end
	z = z(finite(z));	 % Strip infinities from zeros
	pd = (1+p/fs)./(1-p/fs); % Do bilinear transformation
	zd = (1+z/fs)./(1-z/fs);
% real(kd) or just kd?
	kd = (k*prod(fs-z)./prod(fs-p));
	zd = [zd;-ones(length(pd)-length(zd),1)];  % Add extra zeros at -1

elseif (md == 1 & mn == 1) | nargout == 4 %
	if nargout == 4		% State-space case
		a = z; b = p; c = k; d = fs; fs = fp;
		error(abcdchk(a,b,c,d));
		if nargin == 6			% Prewarp
			fp = fp1;		% Decode arguments
			fp = 2*pi*fp;
			fs = fp/tan(fp/fs/2)/2;
		end
	else			% Transfer function case
		if nn > nd
			error('Numerator cannot be higher order than denominator.')
		end
		num = z; den = p;		% Decode arguments
		if nargin == 4			% Prewarp
			fp = fs; fs = k;	% Decode arguments
			fp = 2*pi*fp;
			fs = fp/tan(fp/fs/2)/2;
		else
			fs = k;			% Decode arguments
		end
		% Put num(s)/den(s) in state-space canonical form.  
		[a,b,c,d] = tf2ss(num,den);
	end
	% Now do state-space version of bilinear transformation:
	t = 1/fs;
	r = sqrt(t);
	t1 = eye(size(a)) + a*t/2;
	t2 = eye(size(a)) - a*t/2;
	ad = t2\t1;
	bd = t/r*(t2\b);
	cd = r*c/t2;
	dd = c/t2*b*t/2 + d;
	if nargout == 4
		zd = ad; pd = bd; kd = cd;
	else
		% Convert back to transfer function form:
		p = poly(ad);
		zd = poly(ad-bd*cd)+(dd-1)*p;
		pd = p;
	end
else
	error('First two arguments must have the same orientation.')
end

function [at,bt,ct,dt] = lp2lp(a,b,c,d,wo)
%LP2LP Lowpass to lowpass analog filter transformation.
%   [NUMT,DENT] = LP2LP(NUM,DEN,Wo) transforms the lowpass filter
%   prototype NUM(s)/DEN(s) with unity cutoff frequency to a
%   lowpass filter with cutoff frequency Wo.
%   [AT,BT,CT,DT] = LP2LP(A,B,C,D,Wo) does the same when the
%   filter is described in state-space form.

if nargin == 3		% Transfer function case
        % handle column vector inputs: convert to rows
        if size(a,2) == 1
            a = a(:).';
        end
        if size(b,2) == 1
            b = b(:).';
        end
	% Transform to state-space
	wo = c;
	[a,b,c,d] = tf2ss(a,b);
end

error(abcdchk(a,b,c,d));
[ma,nb] = size(b);
[mc,ma] = size(c);

% Transform lowpass to lowpass
at = wo*a;
bt = wo*b;
ct = c;
dt = d;

if nargin == 3		% Transfer function case
    % Transform back to transfer function
    [z,k] = tzero(at,bt,ct,dt);
    num = k * poly(z);
    den = poly(at);
    at = num;
    bt = den;
end

function [btype,analog,errStr] = iirchk(Wn,varargin)
%IIRCHK  Parameter checking for BUTTER, CHEBY1, CHEBY2, and ELLIP.
%   [btype,analog,errStr] = iirchk(Wn,varargin) returns the 
%   filter type btype (1=lowpass, 2=bandpss, 3=highpass, 4=bandstop)
%   and analog flag analog (0=digital, 1=analog) given the edge
%   frequency Wn (either a one or two element vector) and the
%   optional arguments in varargin.  The variable arguments are 
%   either empty, a one element cell, or a two element cell.
%
%   errStr is empty if no errors are detected; otherwise it contains
%   the error message.  If errStr is not empty, btype and analog
%   are invalid.

errStr = '';

% Define defaults:
analog = 0; % 0=digital, 1=analog
btype = 1;  % 1=lowpass, 2=bandpss, 3=highpass, 4=bandstop

if length(Wn)==1
    btype = 1;  
elseif length(Wn)==2
    btype = 2;
else
    errStr = 'Wn must be a one or two element vector.';
    return
end

if length(varargin)>2
    errStr = 'Too many input arguments.';
    return
end

% Interpret and strip off trailing 's' or 'z' argument:
if length(varargin)>0 
    switch lower(varargin{end})
    case 's'
        analog = 1;
        varargin(end) = [];
    case 'z'
        analog = 0;
        varargin(end) = [];
    otherwise
        if length(varargin) > 1
            errStr = 'Analog flag must be either ''z'' or ''s''.';
            return
        end
    end
end

% At this point, varargin will either be empty, or contain a single
% band type flag.

if length(varargin)==1   % Interpret filter type argument:
    switch lower(varargin{1})
    case 'low'
        btype = 1;
    case 'bandpass'
        btype = 2;
    case 'high'
        btype = 3;
    case 'stop'
        btype = 4;
    otherwise
        if nargin == 2
            errStr = ['Option string must be one of ''high'', ''stop'',' ...
              ' ''low'', ''bandpass'', ''z'' or ''s''.'];
        else  % nargin == 3
            errStr = ['Filter type must be one of ''high'', ''stop'',' ...
              ' ''low'', or ''bandpass''.'];
        end
        return
    end
    switch btype
    case 1
        if length(Wn)~=1
            errStr = 'For the ''low'' filter option, Wn must have 1 element.';
            return
        end
    case 2
        if length(Wn)~=2
            errStr = 'For the ''bandpass'' filter option, Wn must have 2 elements.';
            return
        end
    case 3
        if length(Wn)~=1
            errStr = 'For the ''high'' filter option, Wn must have 1 element.';
            return
        end
    case 4
        if length(Wn)~=2
            errStr = 'For the ''stop'' filter option, Wn must have 2 elements.';
            return
        end
    end
end

function [z,p,k] = besselap(n)
%BESSELAP  Bessel analog lowpass filter prototype.
%   [Z,P,K] = BESSELAP(N) returns the zeros, poles, and gain
%   for an N-th order normalized prototype Bessel analog
%   lowpass filter.  The cutoff or 3dB frequency is equal to 1 
%   for N = 1 and decreases as N increases.
%
%   See also BESSELF, BUTTAP, CHEB1AP, and CHEB2AP.

%   Author(s): T. Krauss, 3-23-93
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.11 $  $Date: 1997/12/02 18:36:20 $

    if (n>25),
        error('Bessel filter roots are not available for order N > 25.')
    end

    z = [];
    k = 1;

% table method
str = '';
if n == 0,
   p = []; return
elseif n == 1,
   str = '-1';
elseif n == 2,
   str      ='[-.8660254037844386467637229+.4999999999999999999999996*i]';
   str( 2,:)='[-.8660254037844386467637229-.4999999999999999999999996*i]';
elseif n == 3,
   str      ='[                             -.9416000265332067855971980]';
   str( 2,:)='[-.7456403858480766441810907-.7113666249728352680992154*i]';
   str( 3,:)='[-.7456403858480766441810907+.7113666249728352680992154*i]';
elseif n == 4,
   str      ='[-.6572111716718829545787781-.8301614350048733772399715*i]';
   str( 2,:)='[-.6572111716718829545787788+.8301614350048733772399715*i]';
   str( 3,:)='[-.9047587967882449459642637-.2709187330038746636700923*i]';
   str( 4,:)='[-.9047587967882449459642624+.2709187330038746636700926*i]';
elseif n == 5,
   str      ='[                             -.9264420773877602247196260]';
   str( 2,:)='[-.8515536193688395541722677-.4427174639443327209850002*i]';
   str( 3,:)='[-.8515536193688395541722677+.4427174639443327209850002*i]';
   str( 4,:)='[-.5905759446119191779319432-.9072067564574549539291747*i]';
   str( 5,:)='[-.5905759446119191779319432+.9072067564574549539291747*i]';
elseif n == 6,
   str      ='[-.9093906830472271808050953-.1856964396793046769246397*i]';
   str( 2,:)='[-.9093906830472271808050953+.1856964396793046769246397*i]';
   str( 3,:)='[-.7996541858328288520243325-.5621717346937317988594118*i]';
   str( 4,:)='[-.7996541858328288520243325+.5621717346937317988594118*i]';
   str( 5,:)='[-.5385526816693109683073792-.9616876881954277199245657*i]';
   str( 6,:)='[-.5385526816693109683073792+.9616876881954277199245657*i]';
elseif n == 7,
   str      ='[                             -.9194871556490290014311619]';
   str( 2,:)='[-.8800029341523374639772340-.3216652762307739398381830*i]';
   str( 3,:)='[-.8800029341523374639772340+.3216652762307739398381830*i]';
   str( 4,:)='[-.7527355434093214462291616-.6504696305522550699212995*i]';
   str( 5,:)='[-.7527355434093214462291616+.6504696305522550699212995*i]';
   str( 6,:)='[-.4966917256672316755024763-1.002508508454420401230220*i]';
   str( 7,:)='[-.4966917256672316755024763+1.002508508454420401230220*i]';
elseif n == 8,
   str      ='[-.9096831546652910216327629-.1412437976671422927888150*i]';
   str( 2,:)='[-.9096831546652910216327629+.1412437976671422927888150*i]';
   str( 3,:)='[-.8473250802359334320103023-.4259017538272934994996429*i]';
   str( 4,:)='[-.8473250802359334320103023+.4259017538272934994996429*i]';
   str( 5,:)='[-.7111381808485399250796172-.7186517314108401705762571*i]';
   str( 6,:)='[-.7111381808485399250796172+.7186517314108401705762571*i]';
   str( 7,:)='[-.4621740412532122027072175-1.034388681126901058116589*i]';
   str( 8,:)='[-.4621740412532122027072175+1.034388681126901058116589*i]';
elseif n == 9,
   str      ='[                             -.9154957797499037686769223]';
   str( 2,:)='[-.8911217017079759323183848-.2526580934582164192308115*i]';
   str( 3,:)='[-.8911217017079759323183848+.2526580934582164192308115*i]';
   str( 4,:)='[-.8148021112269012975514135-.5085815689631499483745341*i]';
   str( 5,:)='[-.8148021112269012975514135+.5085815689631499483745341*i]';
   str( 6,:)='[-.6743622686854761980403401-.7730546212691183706919682*i]';
   str( 7,:)='[-.6743622686854761980403401+.7730546212691183706919682*i]';
   str( 8,:)='[-.4331415561553618854685942-1.060073670135929666774323*i]';
   str( 9,:)='[-.4331415561553618854685942+1.060073670135929666774323*i]';
elseif n == 10,
   str      ='[-.9091347320900502436826431-.1139583137335511169927714*i]';
   str( 2,:)='[-.9091347320900502436826431+.1139583137335511169927714*i]';
   str( 3,:)='[-.8688459641284764527921864-.3430008233766309973110589*i]';
   str( 4,:)='[-.8688459641284764527921864+.3430008233766309973110589*i]';
   str( 5,:)='[-.7837694413101441082655890-.5759147538499947070009852*i]';
   str( 6,:)='[-.7837694413101441082655890+.5759147538499947070009852*i]';
   str( 7,:)='[-.6417513866988316136190854-.8175836167191017226233947*i]';
   str( 8,:)='[-.6417513866988316136190854+.8175836167191017226233947*i]';
   str( 9,:)='[-.4083220732868861566219785-1.081274842819124562037210*i]';
   str(10,:)='[-.4083220732868861566219785+1.081274842819124562037210*i]';
elseif n == 11,
   str      ='[                             -.9129067244518981934637318]';
   str( 2,:)='[-.8963656705721166099815744-.2080480375071031919692341*i]';
   str( 3,:)='[-.8963656705721166099815744+.2080480375071031919692341*i]';
   str( 4,:)='[-.8453044014712962954184557-.4178696917801248292797448*i]';
   str( 5,:)='[-.8453044014712962954184557+.4178696917801248292797448*i]';
   str( 6,:)='[-.7546938934722303128102142-.6319150050721846494520941*i]';
   str( 7,:)='[-.7546938934722303128102142+.6319150050721846494520941*i]';
   str( 8,:)='[-.6126871554915194054182909-.8547813893314764631518509*i]';
   str( 9,:)='[-.6126871554915194054182909+.8547813893314764631518509*i]';
   str(10,:)='[-.3868149510055090879155425-1.099117466763120928733632*i]';
   str(11,:)='[-.3868149510055090879155425+1.099117466763120928733632*i]';
elseif n == 12,
   str      ='[-.9084478234140682638817772-95506365213450398415258360.0E-27*i]';
   str( 2,:)='[-.9084478234140682638817772+95506365213450398415258360.0E-27*i]';
   str( 3,:)='[      -.8802534342016826507901575-.2871779503524226723615457*i]';
   str( 4,:)='[      -.8802534342016826507901575+.2871779503524226723615457*i]';
   str( 5,:)='[      -.8217296939939077285792834-.4810212115100676440620548*i]';
   str( 6,:)='[      -.8217296939939077285792834+.4810212115100676440620548*i]';
   str( 7,:)='[      -.7276681615395159454547013-.6792961178764694160048987*i]';
   str( 8,:)='[      -.7276681615395159454547013+.6792961178764694160048987*i]';
   str( 9,:)='[      -.5866369321861477207528215-.8863772751320727026622149*i]';
   str(10,:)='[      -.5866369321861477207528215+.8863772751320727026622149*i]';
   str(11,:)='[      -.3679640085526312839425808-1.114373575641546257595657*i]';
   str(12,:)='[      -.3679640085526312839425808+1.114373575641546257595657*i]';
elseif n == 13,
   str      ='[                             -.9110914665984182781070663]';
   str( 2,:)='[-.8991314665475196220910718-.1768342956161043620980863*i]';
   str( 3,:)='[-.8991314665475196220910718+.1768342956161043620980863*i]';
   str( 4,:)='[-.8625094198260548711573628-.3547413731172988997754038*i]';
   str( 5,:)='[-.8625094198260548711573628+.3547413731172988997754038*i]';
   str( 6,:)='[-.7987460692470972510394686-.5350752120696801938272504*i]';
   str( 7,:)='[-.7987460692470972510394686+.5350752120696801938272504*i]';
   str( 8,:)='[-.7026234675721275653944062-.7199611890171304131266374*i]';
   str( 9,:)='[-.7026234675721275653944062+.7199611890171304131266374*i]';
   str(10,:)='[-.5631559842430199266325818-.9135900338325109684927731*i]';
   str(11,:)='[-.5631559842430199266325818+.9135900338325109684927731*i]';
   str(12,:)='[-.3512792323389821669401925-1.127591548317705678613239*i]';
   str(13,:)='[-.3512792323389821669401925+1.127591548317705678613239*i]';
elseif n == 14,
   str      ='[-.9077932138396487614720659-82196399419401501888968130.0E-27*i]';
   str( 2,:)='[-.9077932138396487614720659+82196399419401501888968130.0E-27*i]';
   str( 3,:)='[      -.8869506674916445312089167-.2470079178765333183201435*i]';
   str( 4,:)='[      -.8869506674916445312089167+.2470079178765333183201435*i]';
   str( 5,:)='[      -.8441199160909851197897667-.4131653825102692595237260*i]';
   str( 6,:)='[      -.8441199160909851197897667+.4131653825102692595237260*i]';
   str( 7,:)='[      -.7766591387063623897344648-.5819170677377608590492434*i]';
   str( 8,:)='[      -.7766591387063623897344648+.5819170677377608590492434*i]';
   str( 9,:)='[      -.6794256425119233117869491-.7552857305042033418417492*i]';
   str(10,:)='[      -.6794256425119233117869491+.7552857305042033418417492*i]';
   str(11,:)='[      -.5418766775112297376541293-.9373043683516919569183099*i]';
   str(12,:)='[      -.5418766775112297376541293+.9373043683516919569183099*i]';
   str(13,:)='[      -.3363868224902037330610040-1.139172297839859991370924*i]';
   str(14,:)='[      -.3363868224902037330610040+1.139172297839859991370924*i]';
elseif n == 15,
   str      ='[                             -.9097482363849064167228581]';
   str( 2,:)='[-.9006981694176978324932918-.1537681197278439351298882*i]';
   str( 3,:)='[-.9006981694176978324932918+.1537681197278439351298882*i]';
   str( 4,:)='[-.8731264620834984978337843-.3082352470564267657715883*i]';
   str( 5,:)='[-.8731264620834984978337843+.3082352470564267657715883*i]';
   str( 6,:)='[-.8256631452587146506294553-.4642348752734325631275134*i]';
   str( 7,:)='[-.8256631452587146506294553+.4642348752734325631275134*i]';
   str( 8,:)='[-.7556027168970728127850416-.6229396358758267198938604*i]';
   str( 9,:)='[-.7556027168970728127850416+.6229396358758267198938604*i]';
   str(10,:)='[-.6579196593110998676999362-.7862895503722515897065645*i]';
   str(11,:)='[-.6579196593110998676999362+.7862895503722515897065645*i]';
   str(12,:)='[-.5224954069658330616875186-.9581787261092526478889345*i]';
   str(13,:)='[-.5224954069658330616875186+.9581787261092526478889345*i]';
   str(14,:)='[-.3229963059766444287113517-1.149416154583629539665297*i]';
   str(15,:)='[-.3229963059766444287113517+1.149416154583629539665297*i]';
elseif n == 16,
   str      ='[-.9072099595087001356491337-72142113041117326028823950.0E-27*i]';
   str( 2,:)='[-.9072099595087001356491337+72142113041117326028823950.0E-27*i]';
   str( 3,:)='[      -.8911723070323647674780132-.2167089659900576449410059*i]';
   str( 4,:)='[      -.8911723070323647674780132+.2167089659900576449410059*i]';
   str( 5,:)='[      -.8584264231521330481755780-.3621697271802065647661080*i]';
   str( 6,:)='[      -.8584264231521330481755780+.3621697271802065647661080*i]';
   str( 7,:)='[      -.8074790293236003885306146-.5092933751171800179676218*i]';
   str( 8,:)='[      -.8074790293236003885306146+.5092933751171800179676218*i]';
   str( 9,:)='[      -.7356166304713115980927279-.6591950877860393745845254*i]';
   str(10,:)='[      -.7356166304713115980927279+.6591950877860393745845254*i]';
   str(11,:)='[      -.6379502514039066715773828-.8137453537108761895522580*i]';
   str(12,:)='[      -.6379502514039066715773828+.8137453537108761895522580*i]';
   str(13,:)='[      -.5047606444424766743309967-.9767137477799090692947061*i]';
   str(14,:)='[      -.5047606444424766743309967+.9767137477799090692947061*i]';
   str(15,:)='[      -.3108782755645387813283867-1.158552841199330479412225*i]';
   str(16,:)='[      -.3108782755645387813283867+1.158552841199330479412225*i]';
elseif n == 17,
   str      ='[                             -.9087141161336397432860029]';
   str( 2,:)='[-.9016273850787285964692844-.1360267995173024591237303*i]';
   str( 3,:)='[-.9016273850787285964692844+.1360267995173024591237303*i]';
   str( 4,:)='[-.8801100704438627158492165-.2725347156478803885651973*i]';
   str( 5,:)='[-.8801100704438627158492165+.2725347156478803885651973*i]';
   str( 6,:)='[-.8433414495836129204455491-.4100759282910021624185986*i]';
   str( 7,:)='[-.8433414495836129204455491+.4100759282910021624185986*i]';
   str( 8,:)='[-.7897644147799708220288138-.5493724405281088674296232*i]';
   str( 9,:)='[-.7897644147799708220288138+.5493724405281088674296232*i]';
   str(10,:)='[-.7166893842372349049842743-.6914936286393609433305754*i]';
   str(11,:)='[-.7166893842372349049842743+.6914936286393609433305754*i]';
   str(12,:)='[-.6193710717342144521602448-.8382497252826992979368621*i]';
   str(13,:)='[-.6193710717342144521602448+.8382497252826992979368621*i]';
   str(14,:)='[-.4884629337672704194973683-.9932971956316781632345466*i]';
   str(15,:)='[-.4884629337672704194973683+.9932971956316781632345466*i]';
   str(16,:)='[-.2998489459990082015466971-1.166761272925668786676672*i]';
   str(17,:)='[-.2998489459990082015466971+1.166761272925668786676672*i]';
elseif n == 18,
   str      ='[-.9067004324162775554189031-64279241063930693839360680.0E-27*i]';
   str( 2,:)='[-.9067004324162775554189031+64279241063930693839360680.0E-27*i]';
   str( 3,:)='[      -.8939764278132455733032155-.1930374640894758606940586*i]';
   str( 4,:)='[      -.8939764278132455733032155+.1930374640894758606940586*i]';
   str( 5,:)='[      -.8681095503628830078317207-.3224204925163257604931634*i]';
   str( 6,:)='[      -.8681095503628830078317207+.3224204925163257604931634*i]';
   str( 7,:)='[      -.8281885016242836608829018-.4529385697815916950149364*i]';
   str( 8,:)='[      -.8281885016242836608829018+.4529385697815916950149364*i]';
   str( 9,:)='[      -.7726285030739558780127746-.5852778162086640620016316*i]';
   str(10,:)='[      -.7726285030739558780127746+.5852778162086640620016316*i]';
   str(11,:)='[      -.6987821445005273020051878-.7204696509726630531663123*i]';
   str(12,:)='[      -.6987821445005273020051878+.7204696509726630531663123*i]';
   str(13,:)='[      -.6020482668090644386627299-.8602708961893664447167418*i]';
   str(14,:)='[      -.6020482668090644386627299+.8602708961893664447167418*i]';
   str(15,:)='[      -.4734268069916151511140032-1.008234300314801077034158*i]';
   str(16,:)='[      -.4734268069916151511140032+1.008234300314801077034158*i]';
   str(17,:)='[      -.2897592029880489845789953-1.174183010600059128532230*i]';
   str(18,:)='[      -.2897592029880489845789953+1.174183010600059128532230*i]';
elseif n == 19,
   str      ='[                             -.9078934217899404528985092]';
   str( 2,:)='[-.9021937639390660668922536-.1219568381872026517578164*i]';
   str( 3,:)='[-.9021937639390660668922536+.1219568381872026517578164*i]';
   str( 4,:)='[-.8849290585034385274001112-.2442590757549818229026280*i]';
   str( 5,:)='[-.8849290585034385274001112+.2442590757549818229026280*i]';
   str( 6,:)='[-.8555768765618421591093993-.3672925896399872304734923*i]';
   str( 7,:)='[-.8555768765618421591093993+.3672925896399872304734923*i]';
   str( 8,:)='[-.8131725551578197705476160-.4915365035562459055630005*i]';
   str( 9,:)='[-.8131725551578197705476160+.4915365035562459055630005*i]';
   str(10,:)='[-.7561260971541629355231897-.6176483917970178919174173*i]';
   str(11,:)='[-.7561260971541629355231897+.6176483917970178919174173*i]';
   str(12,:)='[-.6818424412912442033411634-.7466272357947761283262338*i]';
   str(13,:)='[-.6818424412912442033411634+.7466272357947761283262338*i]';
   str(14,:)='[-.5858613321217832644813602-.8801817131014566284786759*i]';
   str(15,:)='[-.5858613321217832644813602+.8801817131014566284786759*i]';
   str(16,:)='[-.4595043449730988600785456-1.021768776912671221830298*i]';
   str(17,:)='[-.4595043449730988600785456+1.021768776912671221830298*i]';
   str(18,:)='[-.2804866851439370027628724-1.180931628453291873626003*i]';
   str(19,:)='[-.2804866851439370027628724+1.180931628453291873626003*i]';
elseif n == 20,
   str      ='[-.9062570115576771146523497-57961780277849516990208850.0E-27*i]';
   str( 2,:)='[-.9062570115576771146523497+57961780277849516990208850.0E-27*i]';
   str( 3,:)='[      -.8959150941925768608568248-.1740317175918705058595844*i]';
   str( 4,:)='[      -.8959150941925768608568248+.1740317175918705058595844*i]';
   str( 5,:)='[      -.8749560316673332850673214-.2905559296567908031706902*i]';
   str( 6,:)='[      -.8749560316673332850673214+.2905559296567908031706902*i]';
   str( 7,:)='[      -.8427907479956670633544106-.4078917326291934082132821*i]';
   str( 8,:)='[      -.8427907479956670633544106+.4078917326291934082132821*i]';
   str( 9,:)='[      -.7984251191290606875799876-.5264942388817132427317659*i]';
   str(10,:)='[      -.7984251191290606875799876+.5264942388817132427317659*i]';
   str(11,:)='[      -.7402780309646768991232610-.6469975237605228320268752*i]';
   str(12,:)='[      -.7402780309646768991232610+.6469975237605228320268752*i]';
   str(13,:)='[      -.6658120544829934193890626-.7703721701100763015154510*i]';
   str(14,:)='[      -.6658120544829934193890626+.7703721701100763015154510*i]';
   str(15,:)='[      -.5707026806915714094398061-.8982829066468255593407161*i]';
   str(16,:)='[      -.5707026806915714094398061+.8982829066468255593407161*i]';
   str(17,:)='[      -.4465700698205149555701841-1.034097702560842962315411*i]';
   str(18,:)='[      -.4465700698205149555701841+1.034097702560842962315411*i]';
   str(19,:)='[      -.2719299580251652601727704-1.187099379810885886139638*i]';
   str(20,:)='[      -.2719299580251652601727704+1.187099379810885886139638*i]';
elseif n == 21,
   str      ='[                             -.9072262653142957028884077]';
   str( 2,:)='[-.9025428073192696303995083-.1105252572789856480992275*i]';
   str( 3,:)='[-.9025428073192696303995083+.1105252572789856480992275*i]';
   str( 4,:)='[-.8883808106664449854431605-.2213069215084350419975358*i]';
   str( 5,:)='[-.8883808106664449854431605+.2213069215084350419975358*i]';
   str( 6,:)='[-.8643915813643204553970169-.3326258512522187083009453*i]';
   str( 7,:)='[-.8643915813643204553970169+.3326258512522187083009453*i]';
   str( 8,:)='[-.8299435470674444100273463-.4448177739407956609694059*i]';
   str( 9,:)='[-.8299435470674444100273463+.4448177739407956609694059*i]';
   str(10,:)='[-.7840287980408341576100581-.5583186348022854707564856*i]';
   str(11,:)='[-.7840287980408341576100581+.5583186348022854707564856*i]';
   str(12,:)='[-.7250839687106612822281339-.6737426063024382240549898*i]';
   str(13,:)='[-.7250839687106612822281339+.6737426063024382240549898*i]';
   str(14,:)='[-.6506315378609463397807996-.7920349342629491368548074*i]';
   str(15,:)='[-.6506315378609463397807996+.7920349342629491368548074*i]';
   str(16,:)='[-.5564766488918562465935297-.9148198405846724121600860*i]';
   str(17,:)='[-.5564766488918562465935297+.9148198405846724121600860*i]';
   str(18,:)='[-.4345168906815271799687308-1.045382255856986531461592*i]';
   str(19,:)='[-.4345168906815271799687308+1.045382255856986531461592*i]';
   str(20,:)='[-.2640041595834031147954813-1.192762031948052470183960*i]';
   str(21,:)='[-.2640041595834031147954813+1.192762031948052470183960*i]';
elseif n == 22,
   str      ='[-.9058702269930872551848625-52774908289999045189007100.0E-27*i]';
   str( 2,:)='[-.9058702269930872551848625+52774908289999045189007100.0E-27*i]';
   str( 3,:)='[      -.8972983138153530955952835-.1584351912289865608659759*i]';
   str( 4,:)='[      -.8972983138153530955952835+.1584351912289865608659759*i]';
   str( 5,:)='[      -.8799661455640176154025352-.2644363039201535049656450*i]';
   str( 6,:)='[      -.8799661455640176154025352+.2644363039201535049656450*i]';
   str( 7,:)='[      -.8534754036851687233084587-.3710389319482319823405321*i]';
   str( 8,:)='[      -.8534754036851687233084587+.3710389319482319823405321*i]';
   str( 9,:)='[      -.8171682088462720394344996-.4785619492202780899653575*i]';
   str(10,:)='[      -.8171682088462720394344996+.4785619492202780899653575*i]';
   str(11,:)='[      -.7700332930556816872932937-.5874255426351153211965601*i]';
   str(12,:)='[      -.7700332930556816872932937+.5874255426351153211965601*i]';
   str(13,:)='[      -.7105305456418785989070935-.6982266265924524000098548*i]';
   str(14,:)='[      -.7105305456418785989070935+.6982266265924524000098548*i]';
   str(15,:)='[      -.6362427683267827226840153-.8118875040246347267248508*i]';
   str(16,:)='[      -.6362427683267827226840153+.8118875040246347267248508*i]';
   str(17,:)='[      -.5430983056306302779658129-.9299947824439872998916657*i]';
   str(18,:)='[      -.5430983056306302779658129+.9299947824439872998916657*i]';
   str(19,:)='[      -.4232528745642628461715044-1.055755605227545931204656*i]';
   str(20,:)='[      -.4232528745642628461715044+1.055755605227545931204656*i]';
   str(21,:)='[      -.2566376987939318038016012-1.197982433555213008346532*i]';
   str(22,:)='[      -.2566376987939318038016012+1.197982433555213008346532*i]';
elseif n == 23,
   str      ='[                             -.9066732476324988168207439]';
   str( 2,:)='[-.9027564979912504609412993-.1010534335314045013252480*i]';
   str( 3,:)='[-.9027564979912504609412993+.1010534335314045013252480*i]';
   str( 4,:)='[-.8909283242471251458653994-.2023024699381223418195228*i]';
   str( 5,:)='[-.8909283242471251458653994+.2023024699381223418195228*i]';
   str( 6,:)='[-.8709469395587416239596874-.3039581993950041588888925*i]';
   str( 7,:)='[-.8709469395587416239596874+.3039581993950041588888925*i]';
   str( 8,:)='[-.8423805948021127057054288-.4062657948237602726779246*i]';
   str( 9,:)='[-.8423805948021127057054288+.4062657948237602726779246*i]';
   str(10,:)='[-.8045561642053176205623187-.5095305912227258268309528*i]';
   str(11,:)='[-.8045561642053176205623187+.5095305912227258268309528*i]';
   str(12,:)='[-.7564660146829880581478138-.6141594859476032127216463*i]';
   str(13,:)='[-.7564660146829880581478138+.6141594859476032127216463*i]';
   str(14,:)='[-.6965966033912705387505040-.7207341374753046970247055*i]';
   str(15,:)='[-.6965966033912705387505040+.7207341374753046970247055*i]';
   str(16,:)='[-.6225903228771341778273152-.8301558302812980678845563*i]';
   str(17,:)='[-.6225903228771341778273152+.8301558302812980678845563*i]';
   str(18,:)='[-.5304922463810191698502226-.9439760364018300083750242*i]';
   str(19,:)='[-.5304922463810191698502226+.9439760364018300083750242*i]';
   str(20,:)='[-.4126986617510148836149955-1.065328794475513585531053*i]';
   str(21,:)='[-.4126986617510148836149955+1.065328794475513585531053*i]';
   str(22,:)='[-.2497697202208956030229911-1.202813187870697831365338*i]';
   str(23,:)='[-.2497697202208956030229911+1.202813187870697831365338*i]';
elseif n == 24,
   str      ='[-.9055312363372773709269407-48440066540478700874836350.0E-27*i]';
   str( 2,:)='[-.9055312363372773709269407+48440066540478700874836350.0E-27*i]';
   str( 3,:)='[      -.8983105104397872954053307-.1454056133873610120105857*i]';
   str( 4,:)='[      -.8983105104397872954053307+.1454056133873610120105857*i]';
   str( 5,:)='[      -.8837358034555706623131950-.2426335234401383076544239*i]';
   str( 6,:)='[      -.8837358034555706623131950+.2426335234401383076544239*i]';
   str( 7,:)='[      -.8615278304016353651120610-.3403202112618624773397257*i]';
   str( 8,:)='[      -.8615278304016353651120610+.3403202112618624773397257*i]';
   str( 9,:)='[      -.8312326466813240652679563-.4386985933597305434577492*i]';
   str(10,:)='[      -.8312326466813240652679563+.4386985933597305434577492*i]';
   str(11,:)='[      -.7921695462343492518845446-.5380628490968016700338001*i]';
   str(12,:)='[      -.7921695462343492518845446+.5380628490968016700338001*i]';
   str(13,:)='[      -.7433392285088529449175873-.6388084216222567930378296*i]';
   str(14,:)='[      -.7433392285088529449175873+.6388084216222567930378296*i]';
   str(15,:)='[      -.6832565803536521302816011-.7415032695091650806797753*i]';
   str(16,:)='[      -.6832565803536521302816011+.7415032695091650806797753*i]';
   str(17,:)='[      -.6096221567378335562589532-.8470292433077202380020454*i]';
   str(18,:)='[      -.6096221567378335562589532+.8470292433077202380020454*i]';
   str(19,:)='[      -.5185914574820317343536707-.9569048385259054576937721*i]';
   str(20,:)='[      -.5185914574820317343536707+.9569048385259054576937721*i]';
   str(21,:)='[      -.4027853855197518014786978-1.074195196518674765143729*i]';
   str(22,:)='[      -.4027853855197518014786978+1.074195196518674765143729*i]';
   str(23,:)='[      -.2433481337524869675825448-1.207298683731972524975429*i]';
   str(24,:)='[      -.2433481337524869675825448+1.207298683731972524975429*i]';
elseif n == 25,
   str      ='[                                   -.9062073871811708652496104]';
   str( 2,:)='[-.9028833390228020537142561-93077131185102967450643820.0E-27*i]';
   str( 3,:)='[-.9028833390228020537142561+93077131185102967450643820.0E-27*i]';
   str( 4,:)='[      -.8928551459883548836774529-.1863068969804300712287138*i]';
   str( 5,:)='[      -.8928551459883548836774529+.1863068969804300712287138*i]';
   str( 6,:)='[      -.8759497989677857803656239-.2798521321771408719327250*i]';
   str( 7,:)='[      -.8759497989677857803656239+.2798521321771408719327250*i]';
   str( 8,:)='[      -.8518616886554019782346493-.3738977875907595009446142*i]';
   str( 9,:)='[      -.8518616886554019782346493+.3738977875907595009446142*i]';
   str(10,:)='[      -.8201226043936880253962552-.4686668574656966589020580*i]';
   str(11,:)='[      -.8201226043936880253962552+.4686668574656966589020580*i]';
   str(12,:)='[      -.7800496278186497225905443-.5644441210349710332887354*i]';
   str(13,:)='[      -.7800496278186497225905443+.5644441210349710332887354*i]';
   str(14,:)='[      -.7306549271849967721596735-.6616149647357748681460822*i]';
   str(15,:)='[      -.7306549271849967721596735+.6616149647357748681460822*i]';
   str(16,:)='[      -.6704827128029559528610523-.7607348858167839877987008*i]';
   str(17,:)='[      -.6704827128029559528610523+.7607348858167839877987008*i]';
   str(18,:)='[      -.5972898661335557242320528-.8626676330388028512598538*i]';
   str(19,:)='[      -.5972898661335557242320528+.8626676330388028512598538*i]';
   str(20,:)='[      -.5073362861078468845461362-.9689006305344868494672405*i]';
   str(21,:)='[      -.5073362861078468845461362+.9689006305344868494672405*i]';
   str(22,:)='[      -.3934529878191079606023847-1.082433927173831581956863*i]';
   str(23,:)='[      -.3934529878191079606023847+1.082433927173831581956863*i]';
   str(24,:)='[      -.2373280669322028974199184-1.211476658382565356579418*i]';
   str(25,:)='[      -.2373280669322028974199184+1.211476658382565356579418*i]';
end

for j=1:size(str,1),
   p(j,:) = eval(str(j,:));
end
